Brian E. J. Rose, University at Albany
You really should be looking at The Climate Laboratory book by Brian Rose, where all the same content (and more!) is kept up to date.
Here you are likely to find broken links and broken code.
This document uses the interactive Jupyter notebook format. The notes can be accessed in several different ways:
github at https://github.com/brian-rose/ClimateModeling_coursewareAlso here is a legacy version from 2015.
Many of these notes make use of the climlab package, available at https://github.com/brian-rose/climlab
In [1]:
# Ensure compatibility with Python 2 and 3
from __future__ import print_function, division
This just repeats what we did in the notebook L06_Radiation.ipynb. Refer back there for more details.
In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import climlab
from metpy.plots import SkewT
# Some imports needed to make and display animations
from IPython.display import HTML
from matplotlib import animation
In [3]:
# The NOAA ESRL server is shutdown! January 2019
ncep_url = "http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/"
ncep_air = xr.open_dataset( ncep_url + "pressure/air.mon.1981-2010.ltm.nc",
decode_times=False)
#url = 'http://apdrc.soest.hawaii.edu:80/dods/public_data/Reanalysis_Data/NCEP/NCEP/clima/pressure/air'
#air = xr.open_dataset(url)
## The name of the vertical axis is different than the NOAA ESRL version..
#ncep_air = air.rename({'lev': 'level'})
In [4]:
# Take global, annual average and convert to Kelvin
weight = np.cos(np.deg2rad(ncep_air.lat)) / np.cos(np.deg2rad(ncep_air.lat)).mean(dim='lat')
Tglobal = (ncep_air.air * weight).mean(dim=('lat','lon','time'))
print( Tglobal)
As we did before, we're going to plot these data on a "Skew-T" diagram.
Here we'll make a function to create the diagram, because later we are going to reuse it several times.
In [5]:
def make_skewT():
fig = plt.figure(figsize=(9, 9))
skew = SkewT(fig, rotation=30)
skew.plot(Tglobal.level, Tglobal, color='black', linestyle='-', linewidth=2, label='Observations')
skew.ax.set_ylim(1050, 10)
skew.ax.set_xlim(-90, 45)
# Add the relevant special lines
skew.plot_dry_adiabats(linewidth=0.5)
skew.plot_moist_adiabats(linewidth=0.5)
#skew.plot_mixing_lines()
skew.ax.legend()
skew.ax.set_xlabel('Temperature (degC)', fontsize=14)
skew.ax.set_ylabel('Pressure (hPa)', fontsize=14)
return skew
In [6]:
skew = make_skewT()
We are now going to work on some single-column models of the vertical temperature profile to understand physical factors determining the observed profile.
Models of radiative transfer slice up the atmospheric air column into a series of layer, and calculate the emission and absorption of radiation within each layer.
It's really just a generalization of the model we already looked at:
The concept of radiative equilibrium means that we ignore all methods of heat exchange except for radiation, and ask what temperature profile would exist under that assumption?
We can answer that question by using a radiative transfer model to explicity compute the shortwave and longwave beams, and the warming/cooling of each layer associated with the radiative sources and sinks of energy.
Basically, we reach radiative equilibrium when energy is received and lost through radiation at the same rate in every layer.
Because of the complicated dependence of emission/absorption features on the wavelength of radiation and the different gases, the beam is divided up into many different pieces representing different parts of the electromagnetic spectrum.
We will not look explicitly at this complexity here, but we will use a model that represents these processes at the same level of detail we would in a GCM.
climlab provides two different "GCM-level" radiation codes:
The links above take you to the online climlab documentation.
We're going to use a model called the Rapid Radiative Transfer Model or RRTMG. This is a "serious" and widely-used radiation model, used in many comprehensive GCMs and Numerical Weather Prediction models.
climlab provides an easy-to-use Python wrapper for the RRTMG code.
Before setting up the model, we need some water vapor data.
We're actually going to use the specific humidity field from our CESM control simulation:
In [7]:
datapath = "http://thredds.atmos.albany.edu:8080/thredds/dodsC/cesm/"
atm_control = xr.open_dataset( datapath + "som_1850_f19/som_1850_f19.cam.h0.clim.nc", decode_times=False)
Take global, annual average of the specific humidity:
In [8]:
Qglobal = ((atm_control.Q * atm_control.gw)/atm_control.gw.mean(dim='lat')).mean(dim=('lat','lon','time'))
Qglobal
Out[8]:
In [9]:
fig, ax = plt.subplots()
ax.plot(Qglobal*1000., Qglobal.lev)
ax.invert_yaxis()
ax.set_ylabel('Pressure (hPa)')
ax.set_xlabel('Specific humidity (g/kg)')
ax.grid()
In [10]:
# Make a model on same vertical domain as the GCM
state = climlab.column_state(lev=Qglobal.lev, water_depth=2.5)
state
Out[10]:
In [11]:
radmodel = climlab.radiation.RRTMG(name='Radiation (all gases)',
state=state,
specific_humidity=Qglobal.values,
albedo = 0.25, # this the SURFACE shortwave albedo
timestep = climlab.constants.seconds_per_day,
)
print(radmodel)
Look at a few interesting properties of the model we just created:
In [12]:
# Here's the state dictionary we already created:
radmodel.state
Out[12]:
In [13]:
# Here are the pressure levels in hPa
radmodel.lev
Out[13]:
There is a dictionary called absorber_vmr that holds the volume mixing ratio of all the radiatively active gases in the column"
In [14]:
radmodel.absorber_vmr
Out[14]:
Most are just a single number because they are assumed to be well mixed in the atmosphere.
The exception is ozone, which has a vertical structure taken from observations. Let's plot it
In [15]:
# E.g. the CO2 content (a well-mixed gas) in parts per million
radmodel.absorber_vmr['CO2'] * 1E6
Out[15]:
In [16]:
# here is the data you need for the plot, as a plain numpy arrays:
print(radmodel.lev)
print(radmodel.absorber_vmr['O3'])
In [ ]:
In [ ]:
The other radiatively important gas is of course water vapor, which is stored separately in the specific_humidity attribute:
In [17]:
# specific humidity in kg/kg, on the same pressure axis
radmodel.specific_humidity
Out[17]:
For details you can look at the documentation
In [18]:
for item in radmodel.input:
print(item)
Many of the parameters control the radiative effects of clouds.
But here we should note that the model is initialized with no clouds at all:
In [19]:
radmodel.cldfrac
Out[19]:
Here are the current temperatures (initial condition):
In [20]:
radmodel.Ts
Out[20]:
In [21]:
radmodel.Tatm
Out[21]:
Now let's take a single timestep:
In [22]:
radmodel.step_forward()
In [23]:
radmodel.Ts
Out[23]:
The surface got warmer!
Let's take a look at all the diagnostic information that was generated during that timestep:
Every climlab model has a diagnostics dictionary. Here we are going to check it out as an xarray dataset:
In [24]:
climlab.to_xarray(radmodel.diagnostics)
Out[24]:
The main "job" of a radiative transfer model it to calculate the shortwave and longwave fluxes up and down between each model layer.
For example:
In [25]:
climlab.to_xarray(radmodel.LW_flux_up)
Out[25]:
These are upward longwave fluxes in W/m2.
Why are there 27 data points, when the model has 26 pressure levels?
In [26]:
radmodel.lev
Out[26]:
In [27]:
radmodel.lev_bounds
Out[27]:
The last element of the flux array represents the upward flux from the surface to the first level:
In [28]:
radmodel.LW_flux_up[-1]
Out[28]:
The value is about 390 W m$^{-2}$.
Why?
In [29]:
sigma = 5.67E-8
sigma * 288**4
Out[29]:
The surface temperature was initialized at 288 K, and the surface is treated as very close to a blackbody in the model.
What about the flux from the top layer out to space?
Two ways to access this information:
In [30]:
radmodel.LW_flux_up[0]
Out[30]:
In [31]:
radmodel.OLR
Out[31]:
Of course there is a whole other set of fluxes for the shortwave radiation.
One diagnostic we will often want to look at is the net energy budget at the top of the atmosphere:
In [32]:
radmodel.ASR - radmodel.OLR
Out[32]:
Is the model gaining or losing energy?
In [33]:
while np.abs(radmodel.ASR - radmodel.OLR) > 0.01:
radmodel.step_forward()
Check the energy budget again:
In [34]:
# Check the energy budget again
radmodel.ASR - radmodel.OLR
Out[34]:
Here's a helper function we'll use to add model temperature profiles to our skew-T plot:
In [35]:
def add_profile(skew, model, linestyle='-', color=None):
line = skew.plot(model.lev, model.Tatm - climlab.constants.tempCtoK,
label=model.name, linewidth=2)[0]
skew.plot(1000, model.Ts - climlab.constants.tempCtoK, 'o',
markersize=8, color=line.get_color())
skew.ax.legend()
In [36]:
skew = make_skewT()
add_profile(skew, radmodel)
skew.ax.set_title('Pure radiative equilibrium', fontsize=18);
What do you think about this model -- data comparison?
Models are for experimenting and playing with!
We have just built a single-column radiation model with several different absorbing gases. We can learn about their effects by taking them away.
In [37]:
# Make an exact clone of our existing model
radmodel_noH2O = climlab.process_like(radmodel)
radmodel_noH2O.name = 'Radiation (no H2O)'
print(radmodel_noH2O)
In [38]:
# Here is the water vapor profile we started with
radmodel_noH2O.specific_humidity
Out[38]:
Now get rid of the water entirely!
In [39]:
radmodel_noH2O.specific_humidity *= 0.
In [40]:
radmodel_noH2O.specific_humidity
Out[40]:
In [ ]:
In [41]:
radmodel_noH2O.step_forward()
while np.abs(radmodel_noH2O.ASR - radmodel_noH2O.OLR) > 0.01:
radmodel_noH2O.step_forward()
In [42]:
radmodel_noH2O.ASR - radmodel_noH2O.OLR
Out[42]:
In [43]:
skew = make_skewT()
for model in [radmodel, radmodel_noH2O]:
add_profile(skew, model)
What do you think you can learn from this about the radiative role of water vapor?
In [ ]:
RRTMG radiation model with prescribed profiles of absorbing gases to calculate pure radiative equilibrium temperature profiles.However the really key takeaway message is that none of these radiative equilibrium profiles look much like the observations in the troposphere.
This strongly suggests that other physical processes (aside from radiation) are important in determining the observed temperature profile.
Plotting on the skew-T diagram makes it clear that all the radiative equilibrium profiles are statically unstable near the surface.
The next step is therefore to look at the effects of convective mixing on the temperatures of the surface and lower troposphere.
To make a Radiative-Convective model, we just take a radiation model and couple it to a convection model!
The "convection" model we're going to use here is available as
climlab.convection.ConvectiveAdjustmentIt is a simple process that looks for lapse rates exceeding a critical threshold and performs an instantaneous adjustment that mixes temperatures to the critical lapse rate while conserving energy.
This is a parameterization of the complex, rapid mixing processes that actually occur in an unstable air column!
Here is some code to put this model together in climlab:
In [44]:
# Make a model on same vertical domain as the GCM
state = climlab.column_state(lev=Qglobal.lev, water_depth=2.5)
rad = climlab.radiation.RRTMG(name='Radiation (net)',
state=state,
specific_humidity=Qglobal.values,
timestep = climlab.constants.seconds_per_day,
albedo = 0.25, # tuned to give reasonable ASR for reference cloud-free model
)
conv = climlab.convection.ConvectiveAdjustment(name='Convection',
state=state,
adj_lapse_rate=6.5,
timestep=rad.timestep,)
rcm = climlab.couple([rad, conv], name='Radiative-Convective Model')
In [45]:
print(rcm)
To get some insight into the interaction between radiation and convection, it's useful to look at the adjustment process from a non-equilibrium initial condition.
Let's make an animation!
The code below is complicated but it is mostly for generating the animation. Focus on the results, not the code here.
In [46]:
def get_tendencies(model):
'''Pack all the subprocess tendencies into xarray.Datasets
and convert to units of K / day'''
tendencies_atm = xr.Dataset()
tendencies_sfc = xr.Dataset()
for name, proc, top_proc in climlab.utils.walk.walk_processes(model, topname='Total', topdown=False):
tendencies_atm[name] = proc.tendencies['Tatm'].to_xarray()
tendencies_sfc[name] = proc.tendencies['Ts'].to_xarray()
for tend in [tendencies_atm, tendencies_sfc]:
# convert to K / day
tend *= climlab.constants.seconds_per_day
return tendencies_atm, tendencies_sfc
def initial_figure(model):
fig = plt.figure(figsize=(14,6))
lines = []
skew = SkewT(fig, subplot=(1,2,1), rotation=30)
# plot the observations
skew.plot(Tglobal.level, Tglobal, color='black', linestyle='-', linewidth=2, label='Observations')
lines.append(skew.plot(model.lev, model.Tatm - climlab.constants.tempCtoK,
linestyle='-', linewidth=2, color='C0', label='RC model (all gases)')[0])
skew.ax.legend()
skew.ax.set_ylim(1050, 10)
skew.ax.set_xlim(-60, 75)
# Add the relevant special lines
skew.plot_dry_adiabats(linewidth=0.5)
skew.plot_moist_adiabats(linewidth=0.5)
skew.ax.set_xlabel('Temperature ($^\circ$C)', fontsize=14)
skew.ax.set_ylabel('Pressure (hPa)', fontsize=14)
lines.append(skew.plot(1000, model.Ts - climlab.constants.tempCtoK, 'o',
markersize=8, color='C0', )[0])
ax = fig.add_subplot(1,2,2, sharey=skew.ax)
ax.set_ylim(1050, 10)
ax.set_xlim(-8,8)
ax.grid()
ax.set_xlabel('Temperature tendency ($^\circ$C day$^{-1}$)', fontsize=14)
color_cycle=['g','b','r','y','k']
#color_cycle=['y', 'r', 'b', 'g', 'k']
tendencies_atm, tendencies_sfc = get_tendencies(rcm)
for i, name in enumerate(tendencies_atm.data_vars):
lines.append(ax.plot(tendencies_atm[name], model.lev, label=name, color=color_cycle[i])[0])
for i, name in enumerate(tendencies_sfc.data_vars):
lines.append(ax.plot(tendencies_sfc[name], 1000, 'o', markersize=8, color=color_cycle[i])[0])
ax.legend(loc='center right');
lines.append(skew.ax.text(-100, 50, 'Day {}'.format(int(model.time['days_elapsed'])), fontsize=12))
return fig, lines
def animate(day, model, lines):
lines[0].set_xdata(np.array(model.Tatm)-climlab.constants.tempCtoK)
lines[1].set_xdata(np.array(model.Ts)-climlab.constants.tempCtoK)
#lines[2].set_xdata(np.array(model.q)*1E3)
tendencies_atm, tendencies_sfc = get_tendencies(model)
for i, name in enumerate(tendencies_atm.data_vars):
lines[2+i].set_xdata(tendencies_atm[name])
for i, name in enumerate(tendencies_sfc.data_vars):
lines[2+5+i].set_xdata(tendencies_sfc[name])
lines[-1].set_text('Day {}'.format(int(model.time['days_elapsed'])))
# This is kind of a hack, but without it the initial frame doesn't appear
if day != 0:
model.step_forward()
return lines
We are going to start from an isothermal initial state, and let the model drift toward equilibrium.
In [47]:
# Start from isothermal state
rcm.state.Tatm[:] = rcm.state.Ts
# Call the diagnostics once for initial plotting
rcm.compute_diagnostics()
# Plot initial data
fig, lines = initial_figure(rcm)
Notice several things here:
Now let's look at how the model actually adjusts:
In [48]:
# This is where we make a loop over many timesteps and create an animation in the notebook
ani = animation.FuncAnimation(fig, animate, 150, fargs=(rcm, lines))
HTML(ani.to_html5_video())
Out[48]:
Discuss.
In [49]:
# Here we take JUST THE RADIATION COMPONENT of the full model and run it out to (near) equilibrium
# This is just to get the initial condition for our animation
for n in range(1000):
rcm.subprocess['Radiation (net)'].step_forward()
In [50]:
# Call the diagnostics once for initial plotting
rcm.compute_diagnostics()
# Plot initial data
fig, lines = initial_figure(rcm)
In [51]:
# This is where we make a loop over many timesteps and create an animation in the notebook
ani = animation.FuncAnimation(fig, animate, 100, fargs=(rcm, lines))
HTML(ani.to_html5_video())
Out[51]:
This animation is not as exciting because the instability is destroyed immediately in the first timestep!
That is because the ConvectiveAdjustment process operates instantaneously whenever there is any instability. It is a parameterization taking advantage of the fact that, in nature, convection processes are fast compared to radiative processes.
But notice that the final state is pretty much the same as in the first animation.
The column tends toward the same equilibrium state regardless of where it starts.
Let's repeat our experiment with removing certain absorbing gases from the model, but use the Radiative-Convective model.
In [52]:
# Make a model on same vertical domain as the GCM
state = climlab.column_state(lev=Qglobal.lev, water_depth=2.5)
rad = climlab.radiation.RRTMG(name='all gases',
state=state,
specific_humidity=Qglobal.values,
timestep = climlab.constants.seconds_per_day,
albedo = 0.25, # tuned to give reasonable ASR for reference cloud-free model
)
# remove ozone
rad_noO3 = climlab.process_like(rad)
rad_noO3.absorber_vmr['O3'] *= 0.
rad_noO3.name = 'no O3'
# remove water vapor
rad_noH2O = climlab.process_like(rad)
rad_noH2O.specific_humidity *= 0.
rad_noH2O.name = 'no H2O'
# remove both
rad_noO3_noH2O = climlab.process_like(rad_noO3)
rad_noO3_noH2O.specific_humidity *= 0.
rad_noO3_noH2O.name = 'no O3, no H2O'
# put all models together in a list
rad_models = [rad, rad_noO3, rad_noH2O, rad_noO3_noH2O]
rc_models = []
for r in rad_models:
newrad = climlab.process_like(r)
conv = climlab.convection.ConvectiveAdjustment(name='Convective Adjustment',
state=newrad.state,
adj_lapse_rate=6.5,
timestep=newrad.timestep,)
rc = newrad + conv
rc.name = newrad.name
rc_models.append(rc)
for model in rad_models:
for n in range(100):
model.step_forward()
while (np.abs(model.ASR-model.OLR)>0.01):
model.step_forward()
for model in rc_models:
for n in range(100):
model.step_forward()
while (np.abs(model.ASR-model.OLR)>0.01):
model.step_forward()
In [53]:
skew = make_skewT()
for model in rad_models:
add_profile(skew, model)
skew.ax.set_title('Pure radiative equilibrium', fontsize=18);
In [54]:
skew2 = make_skewT()
for model in rc_models:
add_profile(skew2, model)
skew2.ax.set_title('Radiative-convective equilibrium', fontsize=18);
Lots to discuss here.
The overall message is that equilibrium temperature profile results from a competition between radiation and convection. Essentially:
These calculations all used a critical lapse rate of 6.5 K / km, which is a reasonable approximation to observations.
We set this with the input argument
adj_lapse_rateto the ConvectiveAdjustment process.
The idea is that we are trying to represent the statistical effects of many episodes of moist convection.
An air column that is perfectly neutral to moist instability would follow the blue moist adiabats on the skew-T diagrams.
We can force the model to behave this way by setting
adj_lapse_rate = 'pseudoadiabat'However the real atmosphere, on average, does not exactly follow these adiabats for several reasons:
So here we sweep all this complexity under the rug and just choose a single critical lapse rate for our convective adjustment model.
But this is a parameter that is uncertain and could be interesting to explore.
Repeat the whole series of calculations for the different combinations of absorbing gases above, but with different critical lapse rates:
adj_lapse_rate = 9.8 (in K/km, the dry adiabatic lapse rate, suitable for an atmosphere where condensation does not impact the buoyancy of air parcels)adj_lapse_rate = 'pseudoadiabat' (more suitable for the tropical atmosphere)
In [ ]:
In [55]:
%load_ext version_information
%version_information numpy, matplotlib, xarray, metpy, climlab
Out[55]:
The author of this notebook is Brian E. J. Rose, University at Albany.
It was developed in support of ATM 623: Climate Modeling, a graduate-level course in the Department of Atmospheric and Envionmental Sciences
Development of these notes and the climlab software is partially supported by the National Science Foundation under award AGS-1455071 to Brian Rose. Any opinions, findings, conclusions or recommendations expressed here are mine and do not necessarily reflect the views of the National Science Foundation.
In [ ]: